# R3.6
reference:Controlled modelling of human epiblast and amnion development using stem cells
# loading R library
#R3.6
rm(list=ls())
rewrite=FALSE
condaENV <- "/home/chenzh/miniconda3/envs/R3.6"
LBpath <- paste0(condaENV ,"/lib/R/library")
.libPaths(LBpath)
suppressMessages(library(dplyr))
suppressMessages(library(tidyr))
suppressMessages(library(ggplot2))
suppressMessages(library(cowplot))
suppressMessages(library(pheatmap))
suppressMessages(library(Seurat))
# working directory
DIR <- "~/My_project/CheckBlastoids"
setwd(DIR)
knitr::opts_knit$set(root.dir=DIR)
Loading R functions
#source("~/PC/R_code/functions.R")
#source("~/PC/SnkM/SgCell.R")
source("src/local.quick.fun.R")
feature.plot.cor <- c("yellow","red","black")
suppressMessages(library(foreach))
suppressMessages(library(doParallel))
numCores <- 10
registerDoParallel(numCores)
options(digits = 4)
options(future.globals.maxSize= 3001289600)
TD="Nov14_2021"
loading gene metadata
load("tmp_data/gene.meta.Rdata",verbose=T)
## Loading objects:
## gtf.anno
## mt.gene
## ribo.gene
savefile=paste0("tmp_data/",TD,"/EPSC_Amnion.Rdata")
if (file.exists(savefile)) {
load(savefile,verbose = T)
}else{
load(paste0("tmp_data/",TD,"/all.counts.meta.Rdata"),verbose=T)
#' loading published seurat object
AMLC.ob <- readRDS("data/GSE134571/GSE134571_Posterior48h_H9_Amnion_Merged.rds")
AMLC.ob <- AMLC.ob %>% RenameIdents('0'="Tsw-AMLC",'1'='MeLC2','2'='hES','3'='MeLC1','4'='hPGCLC','5'='AMLC')
AMLC.DM <- FunRF_FindAllMarkers_para(AMLC.ob)
temp.plot1 <- list()
temp.plot1$id <- DimPlot(AMLC.ob,label=T) +NoAxes()+NoLegend()
temp.plot1$oid <- DimPlot(AMLC.ob,label=T,group.by="orig.ident") +NoAxes()+NoLegend()
temp.mk.gene <- (AMLC.DM$sig %>% filter(NP=="pos") %>% filter(gene %in% rownames(counts.all)) %>% group_by(set) %>% top_n(2,power) %>% pull(gene))
for (g in temp.mk.gene) {
temp.plot1[[g]] <- FeaturePlot(AMLC.ob,g)+NoAxes()+NoLegend()+ggtitle(g)+FunTitle()
}
#same QC from the paper
temp.raw <- meta.all %>% filter(pj=="JPF2019")
temp.M <- meta.all %>% filter(pj=="JPF2019" & EML=="Pos48ELS") %>% filter(nGene >=3200 & nGene <= 6200 & mt.perc < 0.06) %>% bind_rows(meta.all %>% filter(pj=="JPF2019" & EML=="H9AMN") %>% filter(nGene >=3600 & nGene <= 6400 & mt.perc < 0.06))
counts.filter <- counts.all[setdiff(rownames(counts.all), mt.gene),temp.M$cell]
DR.TPM.filter <- apply(counts.filter,2,function(x){x/sum(x)*1000000})
# expressed genes
temp.sel.expG <- rownames(counts.filter)[rowSums(counts.filter[,c(temp.M$cell)] >=1) >5]
data.temp <- CreateSeuratObject(counts.filter[temp.sel.expG,temp.M$cell], meta.data = (temp.M %>% tibble::column_to_rownames("cell"))) %>% NormalizeData(verbose = FALSE) %>% FindVariableFeatures( selection.method = "vst", nfeatures = 2000, verbose = FALSE) %>% ScaleData(verbose=F,vars.to.regress=c("mt.perc"))%>% RunPCA(verbose=F,npcs=20) %>% RunUMAP(dims=1:20,verbose=F,min.dist=0.4) %>% FindNeighbors( dims = 1:20,verbose = FALSE) %>% FindClusters(resolution = 0.45,verbose = FALSE) #"mt.perc",
data.mod <- data.temp
id <- Idents(data.mod) %>% as.vector()
names(id) <- colnames(data.mod)
id[colnames(data.mod)[data.mod@meta.data$EML=="H9AMN" & (id %in% c(0,5))]] <- "Tsw-AMLC"
id[colnames(data.mod)[data.mod@meta.data$EML=="H9AMN" & (id %in% c(1))]] <- "hES"
id[colnames(data.mod)[data.mod@meta.data$EML=="Pos48ELS" & (id %in% c(3))]] <- "hPGCLC"
id[colnames(data.mod)[data.mod@meta.data$EML=="Pos48ELS" & (id %in% c(2))]] <- "MeLC1"
id[colnames(data.mod)[data.mod@meta.data$EML=="Pos48ELS" & (id %in% c(4,7))]] <- "MeLC2"
id[colnames(data.mod)[data.mod@meta.data$EML=="Pos48ELS" & (id %in% c(6))]] <- "AMLC"
id[! id %in% c("Tsw-AMLC","hES","hPGCLC","MeLC1","MeLC2","AMLC")] <- "mix"
Idents(data.mod) <- as.factor(id)
data.mod <- subset(data.mod,cells=names(id)[id!="mix"])
table(data.mod %>% Idents())
save(data.mod,temp.sel.expG,temp.M,temp.raw,temp.plot1,AMLC.DM,temp.mk.gene,file=savefile)
}
## Loading objects:
## data.mod
## temp.sel.expG
## temp.M
## temp.raw
## temp.plot1
## AMLC.DM
## temp.mk.gene
cowplot::plot_grid(plotlist=temp.plot1)
temp.plot2 <- list()
temp.plot2$id <- DimPlot(data.mod,label=T) +NoAxes()+NoLegend()
## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.
temp.plot2$oid <- DimPlot(data.mod,label=T,group.by="EML") +NoAxes()+NoLegend()
for (g in temp.mk.gene) {
temp.plot2[[g]] <- FeaturePlot(data.mod,g)+NoAxes()+NoLegend()+ggtitle(g)+FunTitle()
}
cowplot::plot_grid(plotlist=temp.plot2)
#VlnPlot(data.mod,temp.mk.gene)
temp.plot <- list()
temp.plot$AMN1 <- temp.raw %>% filter(EML=="Pos48ELS")%>% ggplot +geom_point(mapping=aes(x=mt.perc,y=nGene,col=EML))+geom_vline(xintercept = 0.06,linetype = "dashed" )+geom_hline(yintercept = 6200,linetype = "dashed")+geom_hline(yintercept = 3200,linetype = "dashed")+theme_classic()+ggtitle("QC for 10X(Pos48ELS)")+FunTitle()
temp.plot$AMN2 <- temp.raw %>% filter(EML=="H9AMN")%>% ggplot +geom_point(mapping=aes(x=mt.perc,y=nGene,col=EML))+geom_vline(xintercept = 0.06,linetype = "dashed" )+geom_hline(yintercept = 6400,linetype = "dashed")+geom_hline(yintercept = 3600,linetype = "dashed")+theme_classic()+ggtitle("QC for 10X(H9AMN)")+FunTitle()
cowplot::plot_grid(plotlist=temp.plot,nrow=2,ncol=2)
AMN.meta <- temp.M %>% filter(cell %in% colnames(data.mod))
AMN.meta$EML <- Idents(data.mod)[AMN.meta$cell] %>% as.vector()
FunMaSF <- function(temp,n) {
set.seed(123)
return(head(temp[sample(nrow(temp)),],n))
}
temp.out <- NULL
meta.AMN.further.ds <- split(AMN.meta,AMN.meta$EML) %>% lapply(function(x){return(FunMaSF(x,100))}) %>% do.call("bind_rows",.)
if (!file.exists(paste0("tmp_data/",TD,"/meta.AMN.further.rds")) | rewrite) {
print("save output")
saveRDS(AMN.meta,paste0("tmp_data/",TD,"/meta.AMN.further.rds"))
saveRDS(meta.AMN.further.ds,file=paste0("tmp_data/",TD,"/meta.AMN.further.ds.rds"))
}
sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server 7.9 (Maipo)
##
## Matrix products: default
## BLAS/LAPACK: /home/chenzh/miniconda3/envs/R3.6/lib/libopenblasp-r0.3.7.so
##
## locale:
## [1] LC_CTYPE=en_US.utf-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.utf-8 LC_COLLATE=en_US.utf-8
## [5] LC_MONETARY=en_US.utf-8 LC_MESSAGES=en_US.utf-8
## [7] LC_PAPER=en_US.utf-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.utf-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] doParallel_1.0.15 iterators_1.0.12 foreach_1.4.8 Seurat_3.1.4
## [5] pheatmap_1.0.12 cowplot_1.0.0 ggplot2_3.2.1 tidyr_1.1.0
## [9] dplyr_1.0.0
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.0-10 Rtsne_0.15 colorspace_1.4-1
## [4] ellipsis_0.3.1 ggridges_0.5.2 farver_2.0.3
## [7] leiden_0.3.3 listenv_0.8.0 npsurv_0.4-0
## [10] ggrepel_0.8.2 mvtnorm_1.0-12 codetools_0.2-16
## [13] splines_3.6.2 mnormt_1.5-6 lsei_1.2-0
## [16] knitr_1.22 TFisher_0.2.0 jsonlite_1.6.1
## [19] ica_1.0-2 cluster_2.0.8 png_0.1-7
## [22] uwot_0.1.5 sctransform_0.2.1 compiler_3.6.2
## [25] httr_1.4.1 Matrix_1.2-17 lazyeval_0.2.2
## [28] htmltools_0.3.6 tools_3.6.2 rsvd_1.0.3
## [31] igraph_1.2.5 gtable_0.3.0 glue_1.4.1
## [34] RANN_2.6.1 reshape2_1.4.4 rappdirs_0.3.1
## [37] Rcpp_1.0.4.6 Biobase_2.46.0 vctrs_0.3.0
## [40] multtest_2.42.0 gdata_2.18.0 ape_5.3
## [43] nlme_3.1-139 gbRd_0.4-11 lmtest_0.9-37
## [46] xfun_0.6 stringr_1.4.0 globals_0.12.5
## [49] lifecycle_0.2.0 irlba_2.3.3 gtools_3.8.2
## [52] future_1.16.0 MASS_7.3-51.3 zoo_1.8-8
## [55] scales_1.1.1 sandwich_2.5-1 RColorBrewer_1.1-2
## [58] yaml_2.2.0 reticulate_1.14 pbapply_1.4-2
## [61] gridExtra_2.3 stringi_1.4.6 highr_0.8
## [64] mutoss_0.1-12 plotrix_3.7-8 caTools_1.18.0
## [67] BiocGenerics_0.32.0 bibtex_0.4.2.2 Rdpack_0.11-1
## [70] rlang_0.4.6 pkgconfig_2.0.3 bitops_1.0-6
## [73] evaluate_0.13 lattice_0.20-38 ROCR_1.0-7
## [76] purrr_0.3.4 labeling_0.3 patchwork_1.0.0
## [79] htmlwidgets_1.3 tidyselect_1.1.0 RcppAnnoy_0.0.15
## [82] plyr_1.8.6 magrittr_1.5 R6_2.4.1
## [85] gplots_3.0.3 generics_0.0.2 multcomp_1.4-12
## [88] pillar_1.4.4 withr_2.2.0 sn_1.5-5
## [91] fitdistrplus_1.0-14 survival_3.1-12 tibble_3.0.1
## [94] future.apply_1.4.0 tsne_0.1-3 crayon_1.3.4
## [97] KernSmooth_2.23-15 plotly_4.9.2 rmarkdown_2.3
## [100] grid_3.6.2 data.table_1.14.0 metap_1.3
## [103] digest_0.6.25 numDeriv_2016.8-1.1 RcppParallel_4.4.4
## [106] stats4_3.6.2 munsell_0.5.0 viridisLite_0.3.0